function result = Gn(k, t, r, s_)
%GN This function gets 4 arguments (k, t, r, s') and returns Gk(t, r)
%   This function implements the following formula:
%   Gk(t, r) = 2pi * integral(-inf => inf)[dt'
%                    integral(0 => inf) [dr' * r'^2
%                    integral(0 => pi) [cos(theta) * Gk-1(t', r') * 
%                                       g(t - t', r, r', theta) d-theta]]]

% Use global parameter:
global u;
global inf_t;
global h_t;

DBG('disp([''k='' num2str(k) ''; t='' num2str(t) ''; r='' num2str(r) ''; s_='' num2str(s_)])');

if (k == 0)
    DBG('disp([''t-(u*s_) = '' num2str(t-(u*s_)) '' len='' num2str(length(t))])');
    result = g(t - (u .* s_), r, 0, 0, s_);
    DBG('disp([''1 Gn res='' num2str(result)])');
else % k >= 1
    %integrand = inline('r_.^2 .* cos(theta) .* Gn(k-1, t_, r_, 500) .* g(t - t_, r, r_, theta, s_)', 'k', 'r', 't', 's_', 't_', 'r_', 'theta');
    %disp(['Gn=' num2str(Gn(k-1, 1, 1, s_))]);
    %disp(['g=' num2str(g(t - 1, r, 1, 1, s_))]);
    %disp(['g*Gn' num2str(Gn(k-1, 1, 1, s_) .* g(t - 1, r, 1, 1, s_))]);
    result = 2 .* pi .* trpzInt(@(t_)int_t(k, r, t, s_, t_), -inf_t, inf_t, h_t);
    %integrand = inline('r_.^2 .* cos(theta) .* Gn(k-1, t_, r_, s_) .* g(t - t_, r, r_, theta, s_)', 'k', 'r', 't', 's_', 't_', 'r_', 'theta');
    %result = 2 .* pi .* triplequad(@(t_, r_, theta)integrand(k, r, t, s_, t_, r_, theta), -10, 10, 0, 10, 0, pi);
    DBG('disp([''2 Gn res='' num2str(result)])');
end

end